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I.  INTENSITY  MODELING 


A.  HISTORICAL  CONTEXT 

The  transmission  of  acoustic  energy  in  waveguides  has  been  a  subject  of  study  for 
many  years.  While  a  complete  review  of  all  progress  to  date  is  beyond  the  scope  of  this 
document,  an  excellent  review  of  the  modeling  and  simulation  techniques  available  to  the 
modern  analyst  is  available  commercially  [1],  While  numerous  scenarios  involve 
radiation  of  energy  into  a  quiescent  free  field,  there  are  similar  problems  in  which  energy 
generated  by  an  acoustic  source  or  scattered  by  a  target  becomes  trapped  by  boundaries 
where  the  description  of  acoustic  propagation  is  within  the  framework  of  guided  waves. 
Examples  of  this  scenario  include  interior  sound  transmission  within  the  air  ventilation 
ducts  of  buildings,  certain  speaker  enclosures  and  most  relevant  to  the  work  discussed 
herein,  cases  of  underwater  transmission.  The  boundaries  acting  to  trap  acoustic  energy 
may  come  from  the  ocean’s  free  surface,  variations  in  the  sound  speed  profile,  and  the 
bathymetry  of  the  environment.  For  example,  sound  is  trapped  in  the  Sound  Fixing  And 
Ranging  (SOFAR)  channel  due  to  fact  that  the  sound  speed  decreases  in  the  thermocline 
to  a  minimum  at  the  SOFAR  axis  and  then  increases  below.  In  this  situation,  the  sound 
energy  is  successively  “bent”  upward  and  downward  while  it  propagates,  not  directly 
interacting  with  a  boundary. 

When  these  are  the  governing  circumstances  of  transmission,  numerical  ray 
tracing  methods  have  historically  been  the  analysis  tool  of  choice.  However, 
these  numerical  techniques  are  sensitive  to  initial  conditions  in  range  dependent 
environments  [2],  and  such  situations  may  necessitate  a  considerable  number  of 
deterministic  trials  to  yield  a  statistically  meaningful  result.  In  cases  where  the  sound 
wave  interacts  with  a  boundary,  there  is  a  predictable,  but  non-trivial  change  in  both  the 
magnitude  and  phase  that  needs  to  be  computed  accurately.  When  this  interaction 
becomes  part  of  the  propagation  path,  ray  methods  may  still  be  applicable,  but  can 
become  cumbersome.  To  address  this,  there  has  been  a  large  body  of  work  devoted  to 
analytical  and  numerical  methods  that  correctly  account  for  the  transmission  of  energy 
within  the  bounded  water  column,  as  well  as  the  interaction  at  the  boundaries.  For 
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brevity,  an  incomplete  listing  of  techniques  includes  methods  of  Normal  Modes  (NM), 
Parabolic  Equations  (PE),  the  aforementioned  ray  tracing  and  the  Finite  Element  Method 
(FEM).  A  summation  of  the  pros  and  cons  will  be  detailed  in  subsequent  sections; 
however,  it  is  worth  noting  that  for  the  environments  studied  herein,  the  two  most 
appropriate  techniques  are  normal  mode  theory  and  the  finite  element  method.  It  is  worth 
noting  that  in  strongly  range  dependent  environments  and  at  higher  frequencies,  the  PE 
method  may  be  more  applicable  and  computationally  efficient.  While  both  techniques  are 
applicable  to  a  wide  range  of  environments,  their  underlying  assumptions  and 
implementation  limitations  make  them  attractive  for  low  frequency,  range -dependent 
problems. 

To  date  there  have  been  several  investigations  into  waveguides  using  normal 
mode  theory,  some  even  supplemented  by  finite  element  techniques,  but  there  has  been 
little  direct  comparison  between  the  two  methods  for  shallow  water,  low  frequency 
propagation.  The  results  of  this  study  compare  the  output  of  a  normal  mode  code,  the 
output  of  a  parabolic  equation  solver  and  the  output  from  the  commercial  finite  element 
software  package  COMSOL  to  examine  the  merits  and  detractors  of  each.  It  is  noted  that 
there  is  generally  excellent  agreement  between  the  finite  element  method  and  normal 
mode  theory,  while  there  are  discrepancies  between  their  results  and  those  predicted  by 
the  parabolic  equation  solver  utilized  here  in  some  environments.  These  discrepancies  are 
attributed  to  the  assumptions  behind  each  method  and  the  specific  implementation  of  the 
boundary  condition  necessary  to  define  the  problem. 

B.  BACKGROUND 

Propagation  of  sound  within  a  waveguide  has  received  much  attention,  dating 
back  to  initial  studies  performed  during  the  Second  World  War.  The  seminal  paper  in  this 
field  was  arguably  authored  by  Pekeris  [3],  who  explored  the  acoustic  pressure  field 
present  within  a  fluid  medium  above  an  infinite  fluid  half  space.  In  the  time  since  there 
have  been  many  additions  to  his  findings.  Although  a  summary  of  additions  to  the  field 
can  be  found  in  the  classic  text  by  Urick  [4],  a  brief  discussion  of  select  recent 
contributions  relevant  to  this  particular  study  follows. 
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Kirby  [5]  studied  the  propagation  of  sound  in  a  waveguide  using  a  hybrid  normal 
mode-finite  element  technique.  His  work  employed  a  mode-matching  scheme  to  couple 
the  detailed  results  of  a  finite  element  calculation  of  the  near  field  generated  by  an 
arbitrarily  complicated  body  to  a  duct  of  constant  cross  section.  Buckingham  [6] 
examined  the  acoustic  field  in  a  Pekeris-type  isospeed  waveguide,  including  the  effects  of 
attenuation  in  the  sediment  layer  over  which  the  acoustic  region  is  located.  His  work  used 
numerically  generated  approximations  for  the  natural  mode  shapes  present  within  the 
water  column.  The  results  indicate  that  higher  order;  dissipative  modes  may  still 
contribute  to  the  measurable  sound  field  when  the  sediment  layer  is  of  sufficiently  high 
sound  speed.  Comparisons  of  normal  mode  techniques  with  experimental  data  taken  in  a 
laboratory  setting  were  carried  out  by  Schneiderwind  [7].  This  work  showed  that  when 
the  shear  and  compressional  elastic  behavior  of  the  sediment  layer  are  significant,  a 
normal  mode  solution  might  be  used  as  an  accurate  predictive  tool  at  ranges  far  from  the 
source. 

The  use  of  finite  element  techniques  has  gained  traction  in  the  past  two  decades 
due  largely  to  the  increase  in  computational  resources  available.  Murphy  [8]  details  the 
formulation  and  implementation  of  the  Finite  Element  Ocean  Model  (FOAM),  which, 
although  unable  to  handle  the  effects  of  shear,  presented  several  advantages  over  the 
normal  mode  and  PE  solutions  at  the  time.  Vendhan  [9]  coupled  the  finite  element 
technique  with  a  normal  mode  based  solution  at  the  terminating  boundaries  and  compared 
his  results  of  a  range  dependent  waveguide  with  those  of  other  numerical  techniques. 

Studies  of  acoustic  intensity,  a  combination  of  both  acoustic  pressure  and  acoustic 
particle  velocity,  have  become  more  common  in  recent  years  with  advances  in  sensor 
technology  and  a  renewed  interest  in  the  potential  merits  of  directional  sensing.  Perhaps 
one  of  the  most  cited  works  in  the  field  of  acoustic  intensity  comes  from  Mann  and  Tichy 
[10].  Their  work  defines  and  explores  the  acoustic  intensity  field  as  not  only  a  real 
component  that  describes  energy  flow,  but  also  an  imaginary  component  that  describes 
energy  storage  within  an  acoustic  medium.  Smith  [11]  subsequently  developed  the 
general  expressions  required  for  computing  the  acoustic  particle  velocity  field  within  the 
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framework  of  normal  mode  and  parabolic  equation  methods  in  range-dependent 
environments. 

More  recently,  Dall’Osto  [12]  examined  the  properties  of  the  intensity  field  for  a 
shallow  water  environment.  This  work  compared  measured  data  to  the  predictions  of  a 
numerical  PE  model  for  a  62m  deep  waveguide.  They  compute  both  the  real  (active)  and 
imaginary  (reactive)  components  of  acoustic  intensity.  They  then  attribute  the  sign  of  the 
active  intensity  to  the  direction  of  propagating  energy  and  are  able  to  identify  the  direct 
and  reflected  contributions  to  the  measured  field.  They  note  that  the  reactive  intensity 
does  not  provide  information  on  the  direction  of  energy  flux,  but  rather  indicate  where 
there  are  spatial  gradients  in  the  acoustic  power  field  (square  magnitude  of  pressure). 
Duan  [13]  presents  a  comparison  of  experimental  data  using  two  different  methods  of 
computing  acoustic  intensity  for  propagation  in  a  waveguide.  The  authors  examine  the 
results  of  a  hybrid  theoretical-numerical  dataset  for  a  non-transient  problem  of 
propagation  in  a  waveguide.  The  method  used  to  generate  theoretical  predictions  is 
complex,  using  a  finite  element  approximation  in  the  near  field  region  immediately 
surrounding  the  source,  to  inform  a  normal  mode  expansion  region  in  the  propagating 
portion  of  the  waveguide  with  eigen  functions  chosen  to  allow  both  forward  and 
backward  propagating  waves  similar  to  the  prior  work  of  Kirby  [5].  Of  note  is  that  the 
experimental  data  to  which  they  make  comparisons  was  performed  in  air  using  a  “p-u” 
type  intensity  probe  (independent,  collocated  pressure  and  velocity  sensors). 

C.  THEORY 

The  three  major  methods  used  to  predict  acoustic  propagation  within  a  shallow 
water  waveguide  at  low  frequencies  are  NM,  PE,  and  FEM.  The  mathematical  foundation 
of  each  method  will  be  briefly  reviewed,  the  coordinate  systems  will  be  defined  and  the 
output  variable  quantities  to  be  compared  will  be  reviewed. 

The  solution  technique  involving  normal  modes  begins  with  the  homogeneous, 
linear,  inviscid  wave  equation  neglecting  any  density  contrasts  at  the  water-sediment 
interface, 
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V2p  +  k2p  =  0, 
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where  p  is  acoustic  pressure  and  k  is  the  acoustic  wavenumber.  In  cylindrical  coordinates, 
a  mathematical  solution  is  sought,  assuming  the  pressure  field  may  be  represented  as  a 
separable  solution  (i.e.,  an  expression  comprised  of  the  product  of  a  functional  form  for 
range  dependence  and  a  functional  form  for  depth  dependence).  For  the  purposes  of 
comparison  between  methods,  it  is  further  assumed  that  there  is  no  range  dependence 
within  the  waveguide.  With  this  simplification,  the  equation  defining  the  normal  modes 
of  pressure  oscillations  within  the  waveguide  is  then  given  as 

d2xVm(z)  ~  2 

— ~j^2 —  =  -FmOO'JmO)  , 

where  represents  the  eigenfunction,  or  mode  shape  for  a  given  modal  index,  m,  and 
ym  is  the  quantized  vertical  component  of  propagating  wavenumber. 

A  solution  to  the  depth-separated  equation  is  then  sought  by  use  of  an  eigenvalue 
solver  in  the  presence  of  certain  boundary  conditions.  For  the  ocean  waveguide,  two 
idealized  boundary  conditions  corresponding  to  a  rigid  bottom  and  pressure  release  free 
surface  may  be  easily  realized.  In  the  case  of  the  pressure  release  surface  (“acoustically 
soft”  boundary),  both  the  acoustic  pressure  and  the  mode  shape  take  on  a  numerical  value 
of  zero,  i.e.. 


vpm(z  =  0)  =  0.  3 

In  the  case  of  an  acoustically  rigid  bottom,  the  normal  component  of  the  acoustic 
particle  velocity  will  be  zero.  This  is  implemented  by  enforcing  the  vertical  derivative  of 
the  mode  shape  to  approach  zero  at  the  depth  of  the  waveguide,  D,  i.e., 

^(z  =  D)  =  0.  4 

This  solution  technique  ensures  that  each  mode  shape  is  unique,  orthogonal  and 

has  real  eigenvalues.  As  will  be  shown  later,  the  ocean  surface  may  be  accurately 

modeled  as  a  true  pressure  release  boundary  at  which  the  pressure  goes  to  zero  when  the 

5 


product  of  the  source  depth  (distance  from  the  free  surface)  and  its  excitation 
wavenumber,  k  are  of  order  0. 1  or  greater.  The  wavenumber,  k,  is  defined  as  the  ratio  of 
the  circular  frequency  (oo)  of  excitation  to  the  sound  speed  (c)  or 

k  ~  5 

c  c 

For  more  general  boundaries,  such  as  those  separating  two  media  of  contrasting 
materials,  the  boundary  condition  becomes  more  complicated.  For  example,  the  boundary 
condition  that  must  be  satisfied  between  a  water  layer,  of  sound  speed  cw  and  density  pw, 
overlying  a  fluid  bottom  layer,  of  sound  speed  q,  and  density  pi,,  takes  the  form 


^m(D)  + 


Pb 


cm 


Pw 


k2_^L_^L 

K  r2  „2 


dz 


■(D)  =  0. 
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The  introduction  of  this  type  of  boundary  condition  creates  several  complex  phenomena, 
as  will  be  described  in  more  detail  later. 

In  contrast  to  the  exact  solution  of  an  analytical  technique  to  the  differential 
equation  describing  the  normal  modes  present  within  a  waveguide,  the  FE  Method  seeks 
an  approximate  but  convergent  solution  to  partial  differential  equations  subject  to 
boundary  conditions  through  the  use  of  linear  algebra  and  polynomial  expansion. 
Solution  to  Equation  (1)  is  sought  by  discretizing  the  problem  domain  into  a  series  of 
smaller  domains  of  finite  size.  These  elements  are  connected  to  each  other  at  nodes,  and 
the  value  of  the  acoustic  pressure  at  these  nodes  is  determined  by  obtaining  the 
coefficient  values  to  a  polynomial  set  of  shape  functions  describing  the  variation  of 
acoustic  pressure  between  nodes.  COMSOL  allows  for  flexibility  in  the  order  of 
polynomial  shape  function  and  permits  the  user  to  choose  first  through  fifth  order. 

Parabolic  Equation  methods  also  begin  with  the  Helmholtz  equation  (Equation  1), 
but  take  a  different  approach  towards  solution.  By  introducing  a  “reduced”  pressure  field, 

1  0 

u(r,z,<f>),  defined  by  p(r,z,</>)  =  —j=u(r,z,</>),  and  defining  the  operator  notation  P  =  — 

Vr  dr 
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the 


and  Q  =(/j  +  e  +  v  +  lf2 ,  where  e  =  n2- 1,  fj  =  -LSL,  and  v  =  1  f  2  , 

k0  3z  k0r  d(j) 

parabolic  equation  formulation  takes  the  form 

(Pop  +  koQop)u  =0.  7 

Using  an  orthogonal  cylindrical  coordinate  system,  the  radial,  azimuthal  and  axial 
ordinates  are  represented  with  variables  r,  z,  and  cp,  respectively.  A  reference  acoustic 
wavenumber  is  defined  as  ko  and  the  medium  specific  index  of  refraction  is  represented  as 
n. 

Assuming  that  acoustic  energy  primarily  propagates  outward  (cylindrically 
divergent  from  the  source  point),  the  Helmholtz  equation  may  be  factored,  leading  to  an 
equation  of  the  parabolic  form  given  by 

.  .  .  .  i  d'¥(r,z,<p)  8 

Pop^Qr,  z,<p)  =  i/c0(?opW,z,<p)  or  - - — - =  Qop'V(r,z>  cp). 

Introducing  a  slowly  modulating  PE  field  function  y/(r,z,</>)  ,  such  thatT  =  i//elk"' ,  the 
pressure  field  may  now  be  represented  by 


p(r,z,(p )  =  P0 


—  i|r(r,  z,  (p)elk°r. 
r 
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This  form  may  be  further  simplified  into  differential  operators  acting  on  the  PE 
field  function,  and  the  treatment  of  these  operators  is  a  field  of  active  research  [14]. 
However,  for  this  course  of  study,  the  Monterey  Miami  Parabolic  Equation  (MMPE) 
formulation  of  Smith  [15],  which  utilizes  a  split-step  Fourier  algorithm,  was  chosen  as 
the  representative  implementation.  Only  the  2-D  version  is  utilized,  since  3-D  effects  are 
beyond  the  scope  of  this  study.  In  addition,  this  form  of  the  MMPE  model  utilized  the 
Thomson-Chapman  wide-angle  approximation  for  the  operator  Qop  [16],  and  employs  the 
traditional  boundary  smoothing  functions  associated  with  the  split-step  Fourier  algorithm. 
The  executable  code  used  is  available  from  the  Ocean  Acoustics  Library  [17]. 
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It  is  worth  noting  that  fundamentally  both  methods  generate  approximate 
solutions  to  the  same  governing  Helmholtz  equation  (Equation  1).  For  normal  mode 
theory  used  herein,  the  assumption  of  range  independence  is  made  explicitly.  Although 
there  are  many  other  implementations  of  normal  mode  theory  (manifest  in  codes,  such  as 
KRAKEN,  COUPLE,  MOATL,  and  others  [17]  that  do  allow  for  range  dependence),  the 
finite  element  method  places  no  such  restrictions  on  the  computational  space.  It  is  also 
worth  noting  that  the  boundary  conditions  to  be  compared  in  the  initial  studies  are  ideal. 
The  choice  of  pressure  release  at  the  free  surface  is  not  mandatory  within  the  normal 
mode  theory,  but  the  relative  impedance  mismatch  between  air  and  water  makes  it 
convenient  to  implement  and  typically  results  in  acceptable  accuracy. 

MMPE  implements  an  image  source  method  to  address  the  free  surface, 
consistent  with  an  ideal  pressure  release  boundary  condition.  While  there  are  potentially 
some  tricks  that  could  be  played  to  simulate  the  air-water  interface  explicitly,  an  ideal 
pressure  release  boundary  is  used  in  all  MMPE  runs.  Similarly,  within  COMSOL  the  free 
surface  may  be  specified  as  an  ideal  pressure  release  boundary,  or  coupled  to  another 
acoustic  media  in  order  to  capture  the  effects  of  reflection  from  a  free  surface  more 
accurately.  Interaction  with  the  bottom  bathymetry  is  an  area  of  active  research,  with 
many  potential  theorems  available  to  address  the  complicated  interaction  between  the 
acoustic  fluid  and  elastic  bottoms  with  varying  degrees  of  stiffness.  Application  of 
boundary  conditions  appropriate  to  reflection  from  sandy,  lossy  sediment  is  possible 
within  both  methods,  but  for  the  purposes  of  comparison,  initial  studies  employ  an  ideal 
acoustically  rigid  bottom.  This  choice  precludes  the  use  of  MMPE  for  initial  studies  due 
to  the  implementation  of  a  blending  function  necessary  to  model  the  water-sediment 
interaction.  For  a  perfectly  rigid  boundary,  the  resulting  impedance  contrast  approaches 
infinity  and  the  MMPE  output  is  entirely  contaminated  by  this  singularity.  However,  in 
later  modeling  efforts  where  the  sediment  is  modeled  as  an  acoustic  media,  all  three 
methods  are  compared  directly. 
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D.  MODELING 


The  simplified  environment  initially  studied  was  an  ideal  two-dimensional 
waveguide  of  constant  bathymetry.  The  water  column  was  chosen  to  be  100  m  deep,  with 
a  50  Hz  continuous  wave  (CW)  omnidirectional  source  placed  at  a  depth  of  75  m  from 
the  free  surface.  A  constant  1500  m/s,  sound  speed  profile  in  a  range  independent  1000 
kg/m3  density  environment  was  used.  In  order  to  examine  the  effects  of  propagation,  the 
waveguide  was  modeled  out  to  a  range  of  500  times  the  depth  of  the  water  column,  or 
5000  m.  In  both  methods,  the  5000  m  range  modeled  was  assumed  to  extend  outward 
towards  infinity,  an  approximation  enforced  in  different  ways  between  the  methods.  In 
the  first  model,  the  bottom  is  assumed  acoustically  rigid,  and  in  the  second,  it  is  a  planar 
interface  with  a  second  acoustic  media  with  sound  speed  of  1900  m/s  and  density  of  1610 
kg/m3.  The  effects  of  shear  are  not  included,  and  there  is  no  attenuation  in  either  acoustic 
media.  A  diagram  of  this  setup  is  shown  in  Figure  1 .  Note  that  p  =  0  represents  a  pressure 
release  surface  and  v  =  0  represents  an  acoustically  rigid  boundary  at  the  bottom  of  the 
waveguide 


P=0 


v=0 

Figure  1 .  Diagram  of  Problem  Setup  (Not  to  Scale) 

The  accuracy  of  any  finite  element  solution  should  be  verified  through  a  mesh 
convergence  study,  but  heuristically  it  has  been  demonstrated  that  8-10  linear  elements 
or  6-8  quadratic  elements  per  wavelength  are  sufficient  to  resolve  propagating  wave 
phenomenon  [18].  While  these  resolution  requirements  are  sufficient  to  recover  the 
pressure  field  accurately,  one  aim  of  this  study  is  also  to  examine  the  intensity  field, 
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which  requires  accurate  calculation  of  the  particle  velocity.  Under  conditions  of  CW 
excitation,  the  particle  velocity  vector  is  related  to  the  harmonically  varying  scalar 
pressure  field  via  spatial  derivatives,  i.e., 

1  dP  10 

Vi  =  - — — . 
loop  OXi 

In  the  finite  element  method,  the  velocity  field  is  determined  by  taking  spatial 
gradients  in  the  direction  of  interest.  The  accurate  recovery  of  gradients  is  strongly  mesh 
dependent,  and  prior  work  has  indicated  that  10-20  quadratic  elements  per  wavelength 
are  required  to  recover  smoothly  varying  spatial  gradients  [19].  Fortunately,  in  the  2D 
axisymmetric  formulation,  for  shallow  water,  low  frequency  scenarios,  this  requirement 
still  allows  for  a  tractable  problem  size. 

In  the  execution  of  the  normal  modes  code,  the  potential  for  examining  a  non- 
rigid  boundary  was  desired.  As  such,  the  formulation  of  Clay  and  Medwin  [20]  was 
coded  into  Matlab  [21].  This  formulation  allows  for  the  interaction  of  the  water  column 
with  a  second  acoustic  media  with  disparate  density  and  sound  speed.  From  a  theoretical 
viewpoint,  the  number  of  modes  present  in  a  waveguide  is  dependent  on  not  only  the 
bounding  geometry,  but  the  boundary  conditions  as  well.  A  thorough  treatment  of  the 
mathematical  techniques  available  for  rigorous  analysis  is  described  by 
Brekhovskikh  [22],  Assuming  cylindrical  spreading,  this  pressure  field  is  the  product  of  a 
radially  propagating  component  knr,  described  by  a  Hankel  function,  and  a  vertical 
standing  component  knz  with  the  total  wavenumber  given  as  the  geometric  mean  of  the 
two,  i.e.. 


k  Jkfir  T  kfiz  . 


A  brief  review  of  the  decomposition  of  the  acoustic  field  present  in  a  waveguide 
with  ideal  boundary  conditions  (rigid  bottom,  pressure  release  top)  can  be  described  as  a 
summation  over  all  modes  that  can  fit  an  odd  number  of  quarter  wavelengths  within  the 
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height,  D,  of  the  waveguide.  Mathematically  this  requires  that 


knz  =  (2n  -  1)  12 

where  n  is  the  modal  index,  and  knz  is  the  axial  component  of  the  acoustic  wavenumber. 

Note  that  because  the  total  wavenumber  is  independent  of  geometry,  there  are 
modes  for  which  the  radial  wavenumber  component  of  propagation  must  become 
imaginary.  When  this  occurs,  the  modes  are  said  to  be  “evanescent.”  Accordingly,  for  a 
given  geometry  and  frequency,  there  are  finite  number  of  modes,  N,  that  will  propagate 
out  radially.  This  can  be  defined  as 

N  =  mod[(^)+i],  13 


Such  modes  are  considered  to  be  “trapped”  because  their  energy  propagates  radially 
outward  from  the  source  at  a  real  angle,  confined  by  the  boundaries  and  are  subject  only 
to  geometric  spreading  and  material/medium  losses,  without  penetration  into  the  other 
bounding  media. 

When  the  water  column  with  sound  speed  c;  sits  atop  a  second  sedimentary 
acoustic  media  with  sound  speed  C2,  the  lower  boundary  condition  requires  that  both  the 
pressure  and  particle  velocities  be  continuous  across  the  interface.  Recall  that  the 
development  of  Transmission  and  Reflection  coefficients  for  incident  plane  waves  at  a 
planar  boundary  require  that  the  incident  and  reflected  angles  as  measured  from  the 
interface  must  be  equal.  Although  there  are  a  number  of  interesting  possibilities,  with 
rigorous  development  and  analysis  available,  [23]  for  the  sake  of  brevity  the  applications 
of  interest  here  are  those  where  the  water  column  has  a  slower  sound  speed.  Under  these 
circumstances  the  critical  angle  ( 0C )  relative  to  the  “grazing  angle”  is  defined  as 


Oc 


-ici 
=  cos  — . 

C2 
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For  all  incident  angles  of  incidence  greater  than  the  critical  angle  there  is  transmitted 

energy  propagating  into  the  sediment  and  the  lower  refraction  angle  is  real,  representing 

1 

the  trapped  modes,  which  decay  in  range  as  ^=.  For  angles  shallower  than  critical,  there  is 

no  transmitted  energy  propagating  into  the  bottom  and  the  lower  refraction  angle 
becomes  complex.  At  the  critical  angle,  there  is  a  singularity  in  the  modal  solution,  which 

represents  a  boundary  wave  in  the  lower  medium  that  propagates  along  the  interface, 
1 

decaying  in  range  as  -.  This  phenomenon  will  lead  to  several  interesting  effects  as  will  be 
shown  later. 

Since  both  media  are  being  excited  at  a  single  frequency,  two  unique  total 

2nf 

wavenumbers  must  exist,  one  in  each  media,  e.g.,  kt  =  — .  The  propagation  vector 

ci 

within  a  given  medium  is  still  dependent  on  the  geometry  as  indicated  by  Equations  (11 
and  12).  Note  however,  the  conditions  of  continuity  also  introduce  the  possibility  of 
transmission  of  energy  across  the  water-sediment  interface,  the  magnitude  of  which 
decreases  exponentially  with  range.  More  precisely,  if  a  mode  corresponds  to  conditions 
involving  propagation  angles  steeper  than  the  critical  angle,  defined  by  Equation  (14),  the 
boundary  condition  itself  becomes  a  complex  valued  function,  and  the  mode,  while  still 
propagating  radially,  decays  exponentially  in  range.  Under  these  circumstances,  energy  is 
propagating  along  the  axis  of  the  waveguide,  but  also  into  the  sediment.  Accordingly, 
energy  in  that  mode  is  no  longer  “trapped”  by  the  boundary  of  the  waveguide,  but  said  to 
be  “leaky.” 

The  number  of  trapped  modes  present  for  a  given  environment  of  depth  D  excited 
at  frequency /with  wave  speeds  in  the  water  column  and  bottom,  c1  and  c2 ,  respectively, 
may  be  calculated  as 


N  =  mod 
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Initial  modeling  attempted  to  mimic  the  effects  of  a  rigid  bottom  by  setting  the  values  of 
sound  speed  and  density  in  the  second  media  to  be  several  orders  of  magnitude  larger 
than  in  the  water  column,  numerically  equal  to  1020  for  both.  In  this  case,  the  second  term 
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in  the  radicand  approaches  zero,  Equation  (15)  approaches  Equation  (13)  and  the  pressure 
field  should  be  represented  by  only  “trapped”  and  not  “leaky”  modes. 

The  pressure  field  excited  by  a  mono  frequency  source  is  then  the  summation  of 
the  product  of  each  trapped  mode  at  all  depths  and  the  value  of  the  mode  at  the  source 
depth,  i.e., 


N 


trapped 


p(r,z)  =  ^ 


771=1 


'Vm(Zs)'VmeiKmr 
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where  m  is  the  modal  index  and  zs  is  the  depth  of  the  source. 

Similarly,  the  velocity  field  components  may  be  obtained  by  taking  spatial 
derivatives  of  the  pressure  field,  and  scaling  by  the  local  impedance  (here  a  constant 
value).  The  complex  intensity  field  may  then  be  simply  computed  by  taking  one-half  the 
product  of  the  pressure  field  and  the  complex  conjugate  of  the  particle  velocity  at  each 
grid  point,  i.e.. 
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where  the  *  denotes  complex  conjugation. 

In  the  simulation  results  presented  here,  a  fixed  grid  with  spacing  of  1  m  in  the 
vertical  direction  and  1  m  in  the  radial  direction  was  used  to  generate  smooth  plots. 
Results  are  not  presented  for  the  sake  of  brevity,  but  through  parametric  variation  of  grid 
spacing,  it  was  found  that  the  results  of  the  normal  mode  code  were  largely  insensitive  to 
output  grid  spacing  for  this  environment. 

Execution  of  the  COMSOL  model  was  carried  out  via  the  pressure  acoustics 
interface  within  version  4.4.  Unlike  the  other  two  methods  where  the  source  may  be 
specified  by  a  delta  function,  the  pressure  field  in  the  absence  of  any  boundaries  may  be 
defined  a  priori  and  the  scattering  interactions  will  be  computed  by  the  software.  This 
method  has  several  advantages,  though  the  implementation  should  be  handled  carefully. 
Within  the  COMSOL  framework,  the  formulation  of  the  Helmholtz  equation  allows  the 
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user  to  solve  for  either  the  “scattered”  or  “total”  pressure  fields.  By  definition,  the 
scattered  field  is  generated  by  interaction  with  boundaries  and/or  scattering  bodies.  The 
incident  field  is  defined  as  the  field,  which  would  exist  in  the  total  absence  of  any 
boundaries  or  scattering  bodies.  The  total  field  is  the  field  that  would  be  measured  by  an 
observer,  and  is  defined  as  the  sum  of  the  scattered  and  incident  fields.  These  definitions 
are  mathematically  viable  since  the  governing  differential  equation  is  linear,  and  each 
field  individually  satisfies  the  governing  equation.  The  sum  of  the  individual  solutions 
will  itself  be  an  admissible  solution.  The  advantages  of  these  definitions  are  twofold: 
first,  it  allows  the  analyst  insight  into  the  relative  position  of  scattered  features  and  the 
source  of  reflections,  and  second,  it  provides  a  natural  separation  of  scales  that  allow  for 
better  conditioned  numerics.  In  general  structural  acoustics  problems  for  which  the  finite 
element  method  is  typically  employed,  the  incident  acoustic  field  may  be  one  or  more 
orders  of  magnitude  larger  than  the  scattered  field.  Since  the  incident  field  is  defined 
before  any  computations  are  carried  out,  only  the  scattered  field  need  be  modeled.  In  the 

case  of  propagation  modeling,  the  field  scattered  from  the  boundaries  is  of  the  same  order 

1 

of  magnitude  as  the  incident  field  (except  immediately  near  the  source  due  to  the  — 

decay),  but  can  be  treated  within  the  same  mathematical  and  numerical  framework.  In 
both  cases,  the  total  field  may  then  be  determined  by  summing  the  scattered  field  with  the 
incident  field.  Implementation  in  COMSOL  is  somewhat  cumbersome,  requiring  that  the 
incident  field  be  defined  mathematically  and  added  to  the  solution  domain  as  a 
“background  pressure  field.” 

Although  the  nomenclature  is  somewhat  misleading,  the  intent  is  to  allow  the 
software  to  compute  accurately  the  relatively  small-scattered  component  of  the  acoustic 
field.  For  the  case  of  the  monopole  source,  this  is  defined  as 

p.  rr  t\  _  I p-i(kr-ut)  18 

where  r  is  the  radial  ordinate,  k  is  the  free  space  acoustic  wavenumber,  and  <x>  is  the 
angular  frequency  of  excitation. 
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Note  that  there  is  a  1/r  decay  in  Equation  (18)  that  results  in  a  singularity  at  the 
acoustic  center  of  the  source.  Finite  element  techniques  cannot  tolerate  such  a  singularity 
within  the  computational  domain,  so  a  scheme  for  dealing  with  this  issue  is  presented  as 
follows.  The  cylindrical  waveguide  is  set-up  in  the  2D  axisymmetric  mode  of  COMSOL. 
In  this  mode,  only  a  single  radial  cross  section  of  the  waveguide  need  be  modeled.  This 
reduces  the  problem  size  down  to  a  tractable  “slice”  of  the  physical  domain,  which  may 
be  recovered  by  revolving  the  computational  domain  about  its  axis.  This,  however,  places 
a  restriction  that  a  point  source  must  be  placed  at  the  axis  of  revolution.  In  order  to 
remove  the  singularity  from  the  computational  domain,  a  circular  portion  of  the  grid, 
centered  at  the  acoustic  center  of  the  point  source  was  removed,  and  an  absorbent 
boundary  placed  on  the  terminating  faces  of  the  resulting  mesh.  COMSOL  has 
implemented  several  methods  for  absorbing  waves,  and  the  use  of  a  “spherical  wave 
radiation”  condition  has  been  shown  to  be  excellent  in  these  situations.  A  closer  view  of 
the  segmented  geometry  with  the  boundary  condition  highlighted  is  shown  in  Figure  2. 


Figure  2.  Close  View  of  Segmented  Geometry  with  Source  Region  and 
Absorbing  Boundary  Highlighted  in  Blue 
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In  the  case  of  the  acoustically  rigid  bottom,  only  the  water  column  must  be 
explicitly  modeled  and  meshed.  At  the  surface  interface,  a  Dirichlet-type  boundary 
condition  is  enforced  whereby  the  pressure  is  equal  to  zero.  The  bottom  interface  can  be 
completely  specified  by  a  Neumann-type  boundary  condition  whereby  the  normal 
derivative  (related  to  the  particle  velocity)  of  the  pressure  field  must  go  to  zero.  Note  that 
because  of  the  relation  between  pressure  gradient  and  particle  velocity,  given  by  equation 
lOthis  is  equivalent  to  enforcing  that  acoustic  particle  motion  match  the  stationary 
boundary  motion,  a  concept  similar  to  the  “no  slip”  boundary  condition  often  employed 
in  the  field  of  fluid  mechanics.  The  complete  geometry  is  represented  by  Figure  1.  In  the 
case  where  the  bottom  is  modeled  as  a  second  acoustic  media,  no  boundary  condition 
need  be  specified  between  the  two  domains;  COMSOL  automatically  enforces  continuity 
of  pressure  and  velocity  along  the  interface.  A  schematic  of  this  situation  is  presented  in 
Figure  3. 
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Figure  3.  Diagram  of  Two  Fluid  Approximation  (Not  to  Scale) 


However,  in  order  to  capture  the  waves  transmitted  across  and  reflected  by  the 
interface,  a  background  pressure  field  must  be  applied  to  the  domain  in  which  the  source 
is  active.  As  stated  above,  the  background  pressure  field  in  the  water  column  is  known  a 
priori  and  may  be  defined  completely  by  a  single  frequency  and  wavenumber.  For  the 

case  of  the  rigid  bottom,  where  only  a  single  sound  speed  is  defined  in  the  water  column, 
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the  choice  of  wavenumber  is  obvious.  In  the  case  where  two  acoustic  media  are  present, 
there  are  two  sound  speeds  and  thus  two  potential  choices  of  wavenumber.  In  the 
COMSOL  framework,  the  interface  between  the  two  domains  acts  as  a  scattering 
boundary.  Because  the  background  field  is  defined  as  that  which  would  be  present  in  the 
absence  of  any  scattering  bodies,  it  is  the  wavenumber  in  the  acoustic  media  in  which  the 
source  is  located  that  must  be  used  to  specify  the  incident  field. 

An  alternate  approach  of  modeling  the  incident  field  was  explored  wherein  a 
spherical  wetted  surface  is  driven  with  constant  radial  acceleration.  This  methodology 
approximates  a  volumetric  monopole  source  and  obviates  the  need  to  define  a 
background  pressure  field,  and  is  consistent  with  historical  theoretical  treatments  of 
monopole  sources.  A  source  wavenumber  must  still  be  specified  and  an  additional  caveat 
that  should  be  noted:  in  order  to  enforce  a  surface  acceleration  such  that  a  constant  source 
level  is  obtained,  the  fluid  properties  and  frequency  of  excitation  must  be  taken  into 
account.  Because  the  acceleration  required  for  a  constant  source  level  is  dependent  on 
both  frequency  and  source  size,  an  equation  based  acceleration,  as  given  in  [24], 


a  (/,  r,R ) 


p0c0  \kr  ) 
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is  enforced  at  all  points  along  the  circular  boundary  defining  the  source  volume. 

The  fluid  properties  (density,  sound  speed  and  wavenumber)  are  again  referenced 
to  the  acoustic  domain  containing  the  source, /is  the  excitation  frequency,  r  is  the  source 
radius  and  p(R)  is  the  pressure  specified  at  some  distance  R.  For  the  studies  presented 
herein,  a  0  deciBel  (dB)  ref  1  pPa  @  lm  source  level  was  desired,  and  the  pressure  at  lm 
from  the  acoustic  center  of  the  finite  size  source  was  held  constant  at  1  pPa.  This 
methodology  was  implemented  as  a  result  of  some  anomalous  findings  and  apparent 
discrepancies  between  the  methods  and  is  discussed  in  detail  in  section  II,  and 
specifically  in  Figure  21.  The  results  are  important  because  they  allow  for  direct 
modeling  of  actual  transduction  devices.  While  the  field  of  fully  coupled  fluid-elastic 
structural  codes  has  been  under  development  for  many  decades,  and  the  demonstration  of 
acceleration  boundaries  to  generate  near  field  results  is  not  new,  the  promise  of  a  fully, 
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directly  coupled  problem  could  yield  new  insight  into  the  governing  mechanisms  of  long 
range  propagation  and  transducer  design  without  the  use  of  intermediary  techniques,  such 
as  the  Kirchoff-Helmholtz  surface  integration  techniques. 

The  incident/scattered  field  separation  approach  has  been  verified  by  modeling  a 
unit  cell  and  enforcing  Floquet-type  [25]  boundary  conditions  between  the  left  and  right 
terminating  faces  to  mimic  an  infinitely  wide  domain.  This  problem  is  created  in  the 
COMSOL  2D  plane  strain  physics  module.  In  this  model,  the  interface  between  the  two 
media  is  modeled  in  the  same  way,  and  both  the  top  and  bottom  boundaries  are 
terminated  with  “planar  wave  radiation”  conditions  to  prevent  spurious  reflections. 
Pointwise  constraints  are  enforced  between  the  left  and  right  boundaries,  separated  by 
distance  L,  such  that  the  horizontal  component  of  the  acoustic  wave-vector  remains 
continuous  across  the  domain  boundaries,  e.g. 

Prim  =  Plefte,k‘L.  20 

A  diagram  showing  the  two  fluids,  location  of  boundaries  and  left  and  right 
pressures  is  shown  in  Figure  4. 


18 


Plane  Wave 
Absorber 


Plane  Wave 
Absorber 

Figure  4.  Diagram  of  Floquet  Periodic  Unit  Cell 

This  cell  is  used  to  verify  transmission  and  reflection  between  two  fluid  domains  using 
the  background  pressure  field  paradigm. 

Here  again,  the  x  component  of  the  wavenumber  used  in  the  constraining  equation 
is  that  of  the  water  column.  A  background  pressure  field  is  defined  in  the  source  acoustic 
medium  by  incident  plane  waves  propagating  in  the  first  fluid  at  some  angle.  In  this  case, 
the  angles  of  incidence  ( 9in )  propagating  in  medium  1  and  refraction  (0re/ra)  in  medium 
2  are  governed  by  Snell’s  law  of  refraction  using  a  coordinate  system  relative  to  normal 
incidence,  mathematically  stated  as 

sin  djji  sin  fly-e/ra  21 

Ci  c2 

This  provides  an  analytic  expression  against  which  the  COMSOL  solution  may  be 
verified.  A  simple  check  on  the  validity  of  this  method  is  performed  using  equal  sound 
speeds  and  densities  in  the  two  acoustic  media  and  observing  that  there  is  no  reflection  or 
refraction  from  the  interface. 
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In  applying  boundary  conditions  to  a  finite  computational  domain  modeling  one 
that  approaches  an  infinitely  large  space,  there  are  several  options  available.  The  first  is 
to  apply  an  “impedance”  boundary  condition.  Experimentation  with  this  option  has 
shown  that  it  prevents  most,  but  not  all  of  the  reflections  at  the  far  end  of  the  domain.  The 
two  more  accurate  options  available  are  the  “radiation”  conditions  and  the  Perfectly 
Matched  Layer  (PML).  The  use  of  PML  has  been  shown  to  be  very  effective  at  higher 
frequencies  [26]  but  requires  an  additional  domain  to  be  meshed  and  several  parameters 
to  be  specified.  In  contrast,  the  radiation  boundary  condition  based  on  the  second  order 
Bayliss-Gunsberger-Turkel  formulation  [27],  designed  to  permit  only  outgoing 
propagating  waves  has  been  shown  to  be  more  effective  at  lower  frequencies. 
Accordingly,  at  the  far  end  of  the  water  column,  a  “cylindrical  wave  radiation”  boundary 
condition  is  applied. 

In  addition  to  the  Pekeris-type  waveguide  [3]  with  a  planar  boundary,  several 
additional  environments  were  studied,  including  deep  water  bathymetries  with  linear 
slopes  converging  to  a  constant  shallow  depth.  This  environment  is  nominally  600  m 
deep,  with  a  section  sloping  upward  to  a  depth  of  150  m  at  a  distance  of  4500  m  from  the 
source,  followed  by  a  second  upward  slope  to  a  depth  of  100  m  at  a  range  of  5500  m  from 
the  source,  finally  ending  in  a  planar  waveguide  of  constant  100  m  depth  for  an  additional 
5000  m  as  shown  in  Figure  5. 


Figure  5.  Schematic  Representation  of  Up  Sloping  Bathymetry 

An  additional  point  of  comparison  between  COMSOL  and  MMPE  comes  from  a 
standard  Pekeris  waveguide  [11].  For  this  study,  the  environment  can  be  summarized  as  a 
200  m  water  column  with  a  constant  sound  speed  of  1500  m/s,  and  density  of  1000  kg/m3 
in  contact  with  a  planar  sediment  layer  possessing  a  constant  sound  speed  of  1700  m/s 
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and  density  of  1500  kg/m3,  excited  by  a  200  Hz  monopole  source  located  at  a  depth  of 
100  m.  Both  fluid  media  are  considered  lossless.  Runtime  parameters  were  provided  by 
the  author  and  execution  of  the  MMPE  code  was  carried  out  locally.  A  comparable 
COMSOL  model  was  created  and  the  transmission  loss  predicted  by  the  two  codes  was 
compared. 

During  the  course  of  modeling  these  environments,  an  additional  point  of 
comparison  and  interesting  comer  case  was  studied,  wherein  the  air-water  interface  was 
not  assumed  a  perfect  reflector  with  known  phase  shift.  Unlike  the  normal  mode 
technique  that  enforces  a  pressure  release  boundary,  or  the  parabolic  equation  technique, 
which  implements  an  image  source  to  model  the  air-water  interface,  COMSOL  can  be 
quickly  configured  to  treat  this  two  fluid  boundary  explicitly,  capturing  all  of  the 
underling  physics.  In  order  to  explore  the  FE  method’s  applicability  in  such 
circumstances,  a  model  of  a  finite  sized  monopole  source  located  in  a  semi-infinite  half 
space  of  air  and  water  contact  was  created.  The  source  was  allowed  to  translate  vertically 
within  the  half-space,  penetrating  the  interface  such  that  the  acoustic  fields  produced  for  a 
continuous  wave  discrete  frequency  excitation  were  available  for  study.  The  boundaries 
of  the  two  half  spaces  were  truncated  with  both  absorbing  boundaries  and  perfectly 
matched  layers  in  order  to  compare  the  efficacy  of  modeling  propagation  from  this  source 
off  to  infinity.  A  schematic  of  the  model  half  spaces  is  presented  in  Figure  6. 
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Figure  6.  Axisymmetric  Model  of  Two  Semi-Infinite  Half  Spaces  with  Planar 

Air-Water  Interface 
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II.  RESULTS 


A.  IDEAL  PEKERIS  WAVEGUIDE 

The  first  set  of  results  comes  from  models  of  the  ideal  rigid  bottom  waveguide. 
Several  relevant  metrics  may  be  used  to  compare  the  output  of  the  two  solution 
methodologies  appropriate  for  this  environment  (finite  element  and  normal  modes).  The 
first  qualitative  method  is  to  examine  color  contour  plots  of  the  complex  acoustic  fields. 
Presented  in  Figure  7  are  contour  plots  of  the  magnitude  of  the  total  acoustic  field  as 
predicted  by  COMSOL  and  the  normal  mode  code. 
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Figure  7.  Comparison  of  COMSOL  and  Normal  Modes:  Acoustic  Magnitude 
Code  predicted  magnitude  of  acoustic  pressure  over  the  0-600  m  range  from  source. 


It  can  be  readily  observed  that  there  is  excellent  qualitative  agreement  between 
the  two  methods.  There  are  some  slight  discrepancies  in  the  location  of  peaks  in  the  field, 
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but  overall  the  features  predicted  by  the  two  methods  are  congruent.  Due  to  the 
formulation  of  the  normal  modes  theory,  and  the  approximations  introduced  by  the 
spherical  absorber  placed  around  the  source  location,  it  is  anticipated  that  the  near  field  of 
both  methods  would  show  minor  differences.  A  comparison  is  therefore  made  to  the 
fields  predicted  by  the  two  codes  at  longer  ranges.  In  addition  to  the  magnitude  of  the 
field,  which  serves  to  illustrate  many  features  of  the  pressure  field,  one  method  to  explore 
the  differences  in  the  two  techniques  is  to  examine  variations  in  complex  phase  angle. 
Therefore,  presented  in  Figure  8  are  contour  plots  of  the  imaginary  component  of  the 
total  pressure  field  over  a  range  of  2000  to  2800  m  from  the  source  location. 
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Figure  8.  Comparison  of  COMSOL  and  Normal  Modes:  Imaginary  Component 

Code  predicted  imaginary  component  of  the  total  acoustic  pressure  field  over  the  2000- 
2800  m  range  from  source. 


At  these  longer  ranges,  the  agreement  between  the  two  codes  is  significantly 
better.  There  are  clearly  identifiable  features  present  in  both  solutions.  As  an  example, 
there  is  a  pronounced  area  of  relatively  high  magnitude  at  2150  m  in  range  and  50  m  in 
depth,  surrounded  by  a  ring  like  set  of  intense  points.  The  surrounding  structure  and 
magnitudes  of  the  features  in  the  pressure  field  has  been  predicted  by  both  methods. 
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While  potentially  interesting  to  examine  in  detail,  the  acoustic  particle  velocity 
field  is  computed  from  spatial  gradients  in  pressure.  Because  of  the  agreement  in  the 
pressure  field,  it  is  anticipated  that  the  velocity  fields  would  agree  as  well.  Rather  than 
illustrate  this  field  directly  and  graphically,  a  comparison  is  made  between  the  two  codes 
of  a  derivative  quantity,  the  acoustic  intensity  (J),  as  defined  by  Equation  17.  For  brevity 
and  a  better  quantitative  metric,  single  line  plots  of  the  active  and  reactive  components  of 
vertical  intensity  at  a  range  of  4500  m  are  presented  in  Figure  9.  The  dotted  red  lines  are 
produced  by  COMSOL,  while  the  continuous  blue  curve  is  produced  by  normal  mode 
theory.  The  leftmost  images  are  the  active  (real)  component  of  vertical  intensity  and  the 
rightmost  images  are  the  reactive  (imaginary)  component. 


x  10-«  Active  Vertical  Intensity  [W/m2]  „  10-«  Reactive  Vertical  Intensity  [W/m2] 


Figure  9.  Line  Plots  of  the  Active  Component  of  the  Vertical  Acoustic  Intensity 
Taken  along  the  depth  of  the  water  column  at  a  range  of  4500  m  from  the  source. 

Examining  the  two  plots,  it  is  clear  that  there  is  excellent  qualitative  and 
quantitative  agreement  between  the  solutions  predicted  by  the  two  methods.  The  shapes 
of  the  vertical  intensity  values  are  consistent,  and  the  magnitudes  predicted  by  both  codes 
are  in  agreement.  There  are  minor  variations  in  the  magnitude  of  the  intensity  peaks 
between  the  two  methods,  but  the  location  of  features  and  their  values  is  consistent. 

It  was  found  that  due  to  the  assumption  of  a  perfectly  reflective  bottom  and  the 
frequencies  studied  in  the  initial  verification  efforts,  the  MMPE  model  was  inappropriate. 

Therefore,  the  results  for  this  environment  are  not  presented.  However,  when 
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environments  that  are  more  realistic  are  modeled,  the  MMPE  technique  produces  high 
quality,  computationally  efficient  predictions  of  the  acoustic  field.  Given  that  the  ultimate 
utility  of  any  model  is  to  predict  the  propagation  of  energy  in  realistic  environments,  and 
all  three  methods  compared  in  this  paper  should,  in  theory,  be  capable  of  handling 
interactions  with  a  second  acoustic  media,  the  results  of  the  environment  described  in  the 
modeling  section  are  now  presented  for  comparison. 

The  verification  of  the  method  of  background  pressure  field  and  acceleration-type 
boundary  is  noteworthy  for  the  two  fluid  domains  because  all  pressure  fields  predicted 
for  a  waveguide  are  predicated  on  the  ability  of  COMSOL  to  handle  the  interface  and 
correctly  predict  several  important  phenomena.  Once  demonstrated  that  the  method  is 
correctly  calculating  the  transmitted  and  reflected  waves  at  this  boundary,  the  results  of  a 
two  fluid  waveguide  are  presented. 

B.  FLOQUET  UNIT  CELL 

Verification  of  the  two  fluid  interaction  is  demonstrated  by  examining  two 
specific  results  from  the  Floquet  unit  cell:  harmonic  plane  waves  incident  on  the  fluid- 
fluid  boundary  at  both  the  sub  critical  angle  and  super  critical  angle.  In  the  case  of 
subcritical  incidence,  the  transmitted  and  reflected  angles  are  both  real  and  readily 
predicted  by  Snell’s  law  as  provided  in  Equation  21.  In  the  case  of  the  supercritical  angle, 
the  transmitted  angle  becomes  imaginary  and  theory  predicts  transition  from  a 
transmitted  wave  to  decaying  evanescent  wave.  In  this  case,  the  amplitude  of  the 
transmitted  wave  decays  exponentially  away  from  the  interface  with  the  energy 
dissipating  primarily  along  the  boundary.  For  the  wave  speeds  chosen  in  this  illustrative 
example,  1500  and  1700  m/s,  the  critical  angle  is  calculated  as  0C  =  62  degrees  relative 
to  vertical.  Presented  in  Figure  10  are  angles  of  incidence  0inc  =  36  degrees  and 
6inc  =  66  degrees.  For  the  sub  critical  incidence  of  36  degrees,  theory  predicts  an  angle 
of  transmission  of  42  degrees. 
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0  me  =  36  [deg]  <  0o  e  inc  =  66  [deg]  >  0c 

Figure  10.  Pressure  Field  Generated  by  Scattering  Off  Fluid  Interface 
At  two  angles  of  incidence,  one  sub  critical  and  one  supercritical. 


In  the  case  of  an  angle  of  incidence  of  36  degrees  shown  on  the  left  of  Figure  10, 
it  can  be  seen  that  energy  propagates  into  the  second  acoustic  media,  and  by  taking  a 
measurement  of  the  angle  of  transmission,  it  can  be  seen  that  the  energy  is  indeed 
propagating  at  42  degrees  relative  to  normal.  In  the  case  of  an  angle  of  incidence  of  66 
degrees  shown  on  the  right  of  Figure  10,  it  can  be  seen  that  the  energy  decays 
exponentially  away  from  the  boundary,  and  does  not  propagate  into  the  second  fluid,  but 
rather  along  the  boundary.  This  also  serves  as  verification  of  the  method  of  implementing 
a  background  pressure  field.  This  plot  is  also  an  excellent  graphical  means  for  illustrating 
evanescent  wave  propagation. 

C.  EXAMINATION  OF  HEAD  WAVE 

While  the  correct  prediction  of  the  evanescent  field  generated  by  interaction  with 
a  second  acoustic  media  has  been  demonstrated,  it  is  necessary  to  verify  the  propagation 
of  the  faster  traveling  wave  in  the  lower  media  to  all  points  in  the  domain.  In  order  to 
show  the  effects  of  the  fast  moving,  laterally  traveling  “head  wave,”  two  visual  aids  are 
presented  below.  First,  a  time  transient  model  was  created  and  driven  at  a  single 
continuous  wave  frequency.  Such  a  model  is  not  well  suited  to  predicting  steady  state 
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phenomenon,  but  does  provide  clear  visual  evidence  of  the  head  wave.  Following  a 
similar  procedure  for  MMPE  execution,  a  schematic  diagram  of  the  water  column, 
sediment  layer  and  deep  bottom  later  is  shown  in  Figure  11.  Note  that  the  boundary 
conditions  between  the  COMSOL  model  shown  in  Figure  3  and  the  MMPE  model  differ 
at  the  terminating  bottom-most  face.  In  the  MMPE  model,  the  deeper  portion  of  the 
sediment  layer  contains  a  spatial  filter  that  acts  to  remove  the  energy  as  it  approaches  the 
edge  of  the  computational  domain,  whereas  in  the  COMSOL  model,  a  plane  wave¬ 
absorbing  boundary  is  used  to  prevent  spurious  reflections  back  to  the  computational 
domain  pressure  field. 
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Figure  11.  Diagram  of  the  Sediment  Bottom  Environment  Modeled  in  MMPE 

(Not  to  Scale) 

The  bathymetry  shown  in  Figure  3  was  used  in  the  aforementioned  time  transient 
run.  The  sound  speed  in  the  second  media  was  parametrically  varied  in  order  to  observe 
the  effects  on  the  acoustic  field  generated  in  the  water  column.  Presented  in  Figure  12  are 
color  contour  plots  of  the  acoustic  fields  predicted  by  COMSOL  for  two  disparate 
sediment  sound  speeds  at  the  same  point  in  simulation  time.  In  the  upper  panel,  the 
sediment  sound  speed  is  set  equal  to  that  of  the  water  column  and  the  early  time,  near 

field  spherical  spreading,  prior  to  any  reflections  from  boundaries,  or  energy  becoming 
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trapped  in  the  cylindrical  section  of  the  waveguide  is  readily  apparent.  Interference 
structure  is  observed  as  reflections  from  the  free  surface  interact  with  the  monopole 
source  field.  In  the  lower  panel,  the  sediment  sound  speed  is  set  to  be  four  times  that  of 
the  water  column.  Under  such  conditions,  the  same  near  field  spherical  spreading  is 
observed.  However,  additional  interference  structure  can  be  seen  in  the  near  field  since 
both  bottom  and  free  surface  reflections  interact  with  the  monopole  source.  More  notably, 
at  farther  range  from  the  source,  the  laterally  traveling  head  wave  can  clearly  be  seen 
moving  away  from  the  source  and  continuously  radiating  acoustic  energy  back  into  the 
water  column  at  an  angle  consistent  with  the  critical  angle  defined  by  Equation  17. 


Case  1  -  c2  =  Ixcl 


Case  2  -  c2  =  4xc1 

Figure  12.  Comparison  of  the  Acoustic  Fields  Predicted  by  COMSOL 
For  two  disparate  sediment  sound  speeds. 


Additional  verification  of  the  correct  prediction  of  the  propagating  head  wave  comes 
from  wavenumber  analysis  of  the  acoustic  field.  The  planar  interface  and  material 
properties  illustrated  in  Figure  3  were  exercised  for  slightly  different  parameters  in  order 
to  provide  adequate  separation  of  the  water  borne  and  sediment  borne  acoustic  waves. 
Specifically,  the  densities  in  both  media  were  set  equal  to  1000  kg/m3  yielding  a  mass 
contrast  ratio,  m  =  —  of  1.0.  The  soundspeed  in  the  water  column  was  set  equal  to  1500 

P2 

m/s  while  the  sediment  was  set  equal  to  2500  m/s,  yielding  an  index  of  refraction,  n  =  — 


of  1.66.  The  monopole  source  was  located  at  a  depth  of  50  m  and  excited  at  a  single 
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frequency  of  100  Hz.  All  surfaces  were  set  to  be  outgoing  wave  absorbers  to  mimic  free- 

space  propagation  of  two  semi-infinite  half  spaces.  Under  such  conditions,  it  is 

anticipated  that  there  should  exist  two  propagating  wave  numbers  with  each 

corresponding  to  the  free -field  motion  of  acoustic  waves.  In  both  the  near  and  far  fields, 

the  water  borne  wave  is  expected  to  propagate  radially  at  /q  =  — .  Similarly,  in  the 

C1 

absence  of  reflection  from  any  boundaries,  the  sediment  wave  is  expected  to  propagate 
radially  at  k2  =  — .  For  the  parameters  listed  above,  these  two  wavenumbers  are 

C2 

numerically  equal  to  /q  =  0.41  rad/m  and  k2  =  0.25  rad/m,  respectively.  By  sequentially 
extracting  short  segments  (sliding  synthetic  aperture  of  4096  points  spaced  at  1  m)  of 
acoustic  data  from  the  pressure  field,  and  spatially  Fourier  transforming  them  into  the 
wavenumber  domain,  via  Matlab’s  FFTW  (Fastest  Fourier  Transform  in  the  West)  [28] 
Discrete  Fourier  Transform  algorithm,  it  is  anticipated  that  two  distinct  peaks  would  be 
observed.  Presented  in  Figure  13  is  a  single  sided  spectrum  obtained  for  one  such  virtual 
aperture. 
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Figure  13.  Single-Sided  Pressure  Amplitude  Spectrum  for 
Sliding  Synthetic  Aperture 

The  Short  Length  Fast  Fourier  Transformed  (SLFFT)  data  is  plotted  as  a 
continuous  blue  curve,  with  green  and  red  circles  indicating  the  theoretical  locations  of 
the  two  expected  propagating  wavenumbers.  Of  note  is  that  the  COMSOL  model 
correctly  predicts  the  pressure  field  and  the  relative  wavenumber  content.  Additionally,  it 
can  be  seen  that  the  head  wave  propagates  at  a  relative  -30  dB  with  respect  to  the  direct 
path.  This  sliding  window  technique  allows  a  sufficiently  long  data  record  to 
accommodate  fine  wavenumber  resolution  from  the  windowed  SLFFT.  By  sampling  the 

acoustic  field  at  dx=l  m  intervals,  the  folding  wavenumber  (similar  to  the  Nyquist  limit 

dx 

of  time  domain  sampling  theory)  is  defined  as  kf0id  =  2n  — ,  which  is  numerically  equal 

to7T  m"1.  Similarly,  the  wavenumber  resolution,  which  corresponds  to  the  bin  spacing  of 

1 

the  SLFFT,  is  given  as  A k  =  2n-^^.  By  choosing  N=4096  (4096  m  for  the  given  dx 

sample  spacing),  this  technique  is  able  to  resolve  wavenumbers  to  within  0.0015  rad/m, 
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which  is  two  orders  of  magnitude  below  the  anticipated  value  of  the  wavenumber 
content. 

By  sequentially  sliding  this  aperture  along  the  length  of  the  water-sediment 
interface,  and  stacking  the  resulting  Fourier  transformed  curves  (as  presented  in  Figure 
13),  a  wavenumber  contour  plot  may  be  constructed  as  a  function  of  the  acoustic  center 
of  the  virtual  aperture.  Performing  this  type  of  processing  allows  for  an  examination  of 
the  propagation  length  scales  and  amplitudes  of  the  acoustic  features  present  for  a  given 
bathymetry.  For  this  simple  case  of  two  semi-infinite  half  spaces,  it  is  anticipated  that 
only  two  dominant  wavenumbers  will  be  present  in  the  far  field.  Figure  14  shows  a 
contour  plot  of  the  wavenumber  content  in  this  environment  as  a  function  of  the  location 
of  the  acoustic  center  of  the  aperture. 
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Figure  14.  Wavenumber  Content  of  the  Two  Semi-Infinite  Half  Spaces 

As  a  function  of  the  location  of  the  acoustic  center  of  the  aperture. 
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It  can  be  seen  that  in  the  very  near  field,  there  is  still  some  spherical  wave  front 
curvature  generated  by  the  monopole  source,  resulting  in  curved  lines  on  the  k-x  contour 
due  to  the  size  of  the  sliding  acoustic  aperture.  At  longer  ranges,  there  are  only  two 
dominant  wavenumbers  propagating,  each  corresponding  with  theoretical  predictions. 
This  type  of  processing  confirms  the  ability  of  COMSOL  to  predict  the  existence  of  the 
evanescent  head  wave. 

The  SLFFT  processing  technique  also  allows  for  a  clear  verification  of  the  effects 
of  sediment  damping  mechanisms.  By  varying  the  value  of  loss  in  the  sediment  from  0.0 
to  1.0  dB/m/kHz,  and  examining  the  propagating  wavenumber  content  in  this  simplified 
geometry,  the  COMSOL  implementation  of  sediment  damping  may  be  illustrated. 
Presented  in  Figure  15  are  three-color  contour  panels  illustrating  the  range  dependence  of 
wavenumber  for  three  unique  damping  values  of  0.0,  0.1,  and  1.0  stacked  from  top  to 
bottom,  respectively. 
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Figure  15.  Range  Dependence  of  Propagating  Wavenumber  Content  within  Semi- 

Infinite  Half  Space 
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In  the  top  panel,  where  the  sediment  damping  is  numerically  equal  to  zero,  the 
lower  valued  wavenumber,  k 2,  corresponding  to  the  sediment  borne  lateral  wave,  is 

consistently  30  dB  less  than  the  direct  path  water  borne  wave,  and  both  wave  types  have 

1 

been  confirmed  to  experience  a  —  spherical  spreading  loss.  As  the  sediment  damping 

value  increases  to  from  0  to  0.1  dB/m/kHz,  it  can  be  seen  in  the  middle  panel  that  the 
head  wave  starts  off  approximately  30  dB  down  relative  to  the  direct  wave.  After 
approximately  3500  m  of  propagation,  the  head  wave  experiences  sufficient  losses  such 
that  it  no  longer  appears  in  the  120  dB  dynamic  range  of  the  figure.  Finally,  in  the  third 
panel,  when  the  sediment  damping  is  increased  to  1.0  dB/m/kHz,  the  head  wave 
contribution  is  damped  out  during  the  initial  1000  m  of  sliding  aperture  such  that  it  no 
longer  appears  in  the  dynamic  range  of  the  plot. 

D.  TWO  FLUID  PEKERIS  WAVEGUIDE 

The  results  of  the  MMPE,  normal  modes,  and  COMSOL  models  are  now 
compared  with  the  inclusion  of  the  lossless  sediment  bottom.  Presented  in  Figure  16  are 
color  contours  of  the  magnitude  of  the  complex  valued  acoustic  field  over  the  range  0- 
600  m  from  the  source.  It  can  be  seen  that  near  the  source,  there  is  reasonable  qualitative 
agreement  between  the  COMSOL  model  and  the  MMPE  model;  however,  there  are 
significant  differences  from  the  normal  modes  model.  In  the  case  of  the  COMSOL  and 
MMPE  fields,  there  appear  to  be  four  “finger”  like  structures  present  in  the  acoustic  field, 
whereas  the  normal  modes  model  does  not  show  any  particular  structure.  The  reason  for 
this  discrepancy  lies  in  the  consideration  of  modes  used  in  the  summation,  and  will  be 
discussed  later  in  the  analysis  section  of  this  document. 
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Figure  16.  Comparison  of  Normal  Modes,  COMSOL  and  MMPE 

Code  predicted  complex  magnitude  of  acoustic  pressure  over  the  0-600  m  range  from 
source  with  the  inclusion  of  a  lossless  sediment  bottom. 

Presentation  of  the  pressure  field  may  reveal  additional  detail  and  discrepancy  if 
the  complex  field  is  broken  into  its  real  and  imaginary  components  and  examined 
separately.  At  farther  ranges,  the  comparison  between  the  codes  improves.  Shown  in 
Figure  17  are  contours  of  the  real  component  of  the  acoustic  pressure  field  over  the 
1000-1600  m  range  from  the  source. 
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Figure  17.  Comparison  of  Normal  Modes,  COMSOL  and  MMPE 

Code  predicted  real  component  of  acoustic  pressure  over  the  0-600  m  range  from  source 
with  the  inclusion  of  a  lossless  sediment  bottom. 


At  these  far  field  ranges,  there  is  a  clear  pattern  in  the  real  component  of  pressure 
that  is  common  to  all  three  codes,  indicative  of  the  time-averaged  interaction  between  the 
limited  number  of  propagating  wave  fronts.  It  is  apparent  that  the  levels  between  the 
three  methods  differ  slightly,  with  the  normal  modes  solution  having  generally  higher 
predicted  values  of  acoustic  pressure  than  the  other  two  codes.  However,  there  is  a 
banded  structure  dominant  in  radial  variation  with  local  minima  that  vary  in  range  over 
the  1100-1150  m  range.  Turning  towards  the  acoustic  intensity  predicted  by  the  three 
codes,  there  are  similar  trends.  Presented  in  Figure  18  are  color  contour  plots  of  the  real 
component  of  the  axial  acoustic  intensity  as  predicted  by  the  three  methods. 


37 


Normal  Modes 


-0.2 


COMSOL 


MMPE 


Range  [m] 

Figure  18.  Contours  of  the  Real  Component  of  the  Axial  Intensity 
Predicted  by  the  three  methods  over  a  1000-1600  m  range  from  the  source. 


Here  again,  the  normal  modes  method  predicts  slightly  higher  levels  than  the 
other  two  codes.  However,  the  general  location  of  regions  of  positive  and  negative 
intensity  is  in  good  qualitative  agreement;  for  instance,  there  is  a  change  of  sign  in  axial 
intensity  predicted  by  all  methods  in  the  1100  to  1300  m  distance  from  the  source.  Also 
consistent  with  other  metrics  of  comparison  is  that  COMSOL  predicts  significantly  more 
structure  in  the  intensity  field  than  do  the  other  two  methods.  The  reasons  will  be 
examined  in  greater  detail  later,  but  this  is  due  to  the  inclusion  of  effects  of  large  angles 
of  propagation. 
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E.  UPSLOPING  BATHYMETRY 


Results  from  the  up  sloping  environment  brought  out  several  differences  between 
solution  techniques.  Discussion  with  NPS  faculty  led  to  the  selection  of  an  academically 
developed  normal  modes  code,  Couple07  [29],  surmised  to  be  capable  of  accurately 
modeling  environments  of  this  type  as  an  alternate  solver. 

Couple07  is  a  2007  update  to  the  FORTRAN  research  code  developed  by  Dr. 
Richard  Evans  that  performs  a  normal  mode  based  calculation  of  the  pressure  field.  This 
code  differs  from  others  in  that  it  allows  for  both  forward  and  backward  propagation 
(coupling)  of  modes  induced  by  interactions  with  downrange  bathymetric  changes  and 
range  dependent  sound  velocity  profiles. 

Therefore,  results  from  Couple  07,  MMPE  and  COMSOL  are  presented  to 
illustrate  the  differences  in  these  techniques.  Figure  19  shows  contour  plots  of  the  active 
vertical  intensity  in  the  up  sloping  environment  for  a  monopole  source  pulsating  at  a 
frequency  of  20  Hz  and  a  depth  of  5  m. 
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Figure  19.  Comparison  of  MMPE,  Couple  07  and  COMSOL 
In  up  sloping  environment  with  a  20  Hz,  5  m  deep  monopole  source. 


Superimposed  on  the  MMPE  results  are  predicted  streamlines  indicating  the  path 
of  energy  flow  in  this  environment.  It  can  be  seen  that  there  is  general  agreement  between 
the  methods  for  this  combination  of  source  depth  and  frequency.  There  is  measurable 
penetration  of  acoustic  intensity  into  the  sediment  at  a  range  of  approximately  3500  m, 
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borne  out  by  the  clustering  of  streamlines.  Increasing  the  frequency  of  excitation  to  60 
Hz,  the  field  produced  by  COMSOL  when  using  the  background  pressure  field 
implementation  shows  a  noticeable  discrepancy  from  the  other  two  techniques  as  shown 
in  Figure  20.  This  is  an  anomalous  prediction  insofar  as  general  agreement  was  found 
when  representing  the  source  via  the  specified  background  field,  volumetric,  or 
acceleration  boundary  condition. 
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Figure  20.  Comparison  of  MMPE,  Couple  07  and  COMSOL 
In  up  sloping  environment  with  a  60  Hz,  5  m  deep  monopole  source. 
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When  the  source  is  represented  instead  with  the  acceleration-type  boundary 
described  earlier,  the  results  from  COMSOL  are  in  much  better  agreement  with  the  other 
two  models.  A  contour  plot  of  the  active  vertical  intensity  field  predicted  by  COMSOL 
using  both  types  of  source  models  is  presented  in  Figure  21. 
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Figure  21.  Comparison  of  the  Source  Model  Types  in  COMSOL 
For  a  60  Hz  monopole  source  operating  at  a  depth  of  5  m. 


After  discovering  the  sensitivity  to  the  source  type  for  this  particular  environment, 
it  was  determined  that  the  acceleration  source  was  the  correct  implementation  of  a 
monopole  under  such  conditions  due  to  the  agreement  with  the  other  two  methodologies 
and  was  thus  used  for  the  remaining  modeling  efforts.  With  this  implementation  of  the 
source,  several  other  frequency-depth  combinations  were  studied.  Presented  in  Figure  22 
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are  the  predictions  for  the  up  sloping  bathymetry  for  a  60  Hz  source  operating  at  a  depth 
of  200  m. 
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Figure  22.  Comparison  of  MMPE,  Couple  07  and  COMSOL  Active  Vertical 

Intensities 

In  up  sloping  environment  with  a  60Hz,  200  m  deep  monopole  source. 
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Again,  there  is  excellent  agreement  between  the  three  codes  in  both  the  location 
of  regions  of  positive  and  negative  intensity,  and  overall  magnitude  of  the  fields.  Due  to 
slight  discrepancies  in  the  scaling  parameters  used,  it  is  more  instructive  to  examine  line 
plots  of  the  transmission  loss  (TL)  defined  as 

TL=20(ofllo(^),  22 

where  pref  =  1  pPa.  Presented  in  Figure  23  are  predictions  of  the  transmission  loss  from 
each  of  the  three  codes  as  a  function  of  range  measured  at  50  m  from  the  free  surface  for 
a  60  Hz  monopole  source  located  at  a  depth  of  200  m. 
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Figure  23.  Transmission  Loss  Plotted  as  a  Function  of  Range 
In  the  up  sloping  environment  for  a  60  Hz  monopole  source  located  at  a  depth  of  200  m. 


As  shown,  there  is  excellent  agreement  between  the  three  codes,  with  the  largest 
discrepancy  being  a  notable  deviation  in  the  phase  predicted  by  MMPE,  a  well- 
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documented  phenomenon,  most  recently  described  by  Erdim  [30],  The  phase  error  is 
generated  during  successive  forward  and  inverse  Fourier  transforms  in  the  presence  of  a 
discontinuity  between  the  two  acoustic  media.  In  order  to  prevent  severe  numerical 
difficulties  several  options  are  available,  with  the  MMPE  code  making  use  of  a  density 
blending  function.  A  slight  (~  2  dB)  difference  between  the  COMSOL  solution  and  the 
other  two  techniques  is  attributed  to  the  accumulated  phase  discrepancies. 

Finally,  comparisons  to  the  results  from  the  Smith  JCA  paper  [11]  are  presented. 
Shown  in  Figure  24  is  a  plot  of  the  transmission  loss  predicted  by  both  COMSOL  and 
MMPE. 


TL  at77.5391m  Depth 


Figure  24.  Comparison  of  the  Predicted  Transmission  Loss  from  COMSOL  and 

MMPE  for  the  Environment  Specified  in  JCA 

It  can  be  seen  that  the  overall  agreement  between  the  two  codes  is  quite  good. 
Consistent  with  prior  comparisons  in  the  shallower  environment,  COMSOL  predicts 
significantly  more  structure  in  the  pressure  field,  especially  at  short  ranges  from  the 
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source.  However,  the  location  of  local  minima  and  maxima  are  in  good  agreement,  as  are 
the  overall  levels  with  increasing  range. 


F.  ANOMALOUS  SURFACE  TRANSPARENCY 

Another  interesting  feature  that  can  be  used  to  investigate  the  accuracy  of  the 
COMSOL  model  was  proposed  by  a  third  party  reviewer  and  relates  to  the  surface 
reflection  from  a  true  air-water  interface  (rather  than  an  ideal  pressure  release 
boundary)  [31].  In  order  to  examine  this  effect,  a  model  was  constructed  with  explicit 
treatment  of  the  two  fluid  media  with  both  sound  speed  and  density  contrasts,  n  and  m, 
respectively.  By  varying  the  vertical  location  of  a  finite  sized  source  within  the  two 
media,  several  interesting  phenomenon  are  observed. 

First,  it  is  noted  that  when  the  source  is  deep  compared  to  an  acoustic  wavelength, 
the  interface  behaves  according  to  classical  theory  (i.e.,  a  nearly  perfect  reflector  of 
acoustic  energy).  Specifically,  when  the  non-dimensional  source  depth  within  the  water 
half  space,  kz,  is  of  order  1  or  greater,  most  of  the  acoustic  energy  remains  in  the  water 
domain.  However,  when  kz  drops  to  order  0.1  or  less,  most  of  the  acoustic  energy  crosses 
the  air-water  boundary  and  becomes  airborne.  Presented  in  several  panels  of  Figure  25 
are  color  contours  of  the  active  vertical  acoustic  intensity  for  positive  and  negative  (water 
borne  and  airborne)  non-dimensional  source  depths  of  order  0.1  and  1. 
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Figure  25.  Contour  Plots  of  Active  Vertical  Intensity  for  Several  Source  Depths 


From  the  color  contours  it  can  be  seen  that  not  only  does  the  majority  of  the 
acoustic  energy  remain  in  the  water  domain  for  large  positive  non-dimensional  source 
depths,  but  what  little  energy  is  transmitted  into  the  air  remains  confined  to  a  narrow  cone 
directly  above  the  source.  At  very  shallow  non-dimensional  source  depths  (order  0.1  or 
less),  there  is  a  sudden  increase  in  the  energy  transmission,  and  the  radiation  pattern 
becomes  increasingly  omnidirectional.  The  relative  fraction  of  acoustic  intensity 
transmitted  across  the  boundary  to  that  reflected  back  into  the  water  domain  is  best 
characterized  by  a  quantity  defined  by  Godin  as  the  Transparency,  T  [31].  Consistent 
with  his  findings,  T  exhibits  a  fundamental  asymmetry  with  respect  to  the  source  depth. 
Presented  in  Figure  26  is  this  ratio  as  a  function  of  non-dimensional  source  depth  as 
predicted  by  the  theory  of  Godin  (left)  and  COMSOL  model  results  (right). 
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Figure  26.  Acoustic  Transparency  of  the  Air-Water  Interface 
As  predicted  by  Godin  (left)  and  COMSOL  (right). 


While  the  predicted  curves  of  Godin  have  been  generated  for  several  values  of  n 
and  m,  corresponding  to  the  curves  labeled  1-3  on  the  left  plot,  the  blue  curve 
corresponds  to  values  typical  of  air  and  water.  For  these  parameter  values,  it  is  noted  that 
there  is  approximately  a  20  dB  difference  between  positive  and  negative  source  depths  of 
order  unity,  while  there  is  nearly  a  50  dB  difference  for  non-dimensional  source  depths  of 
order  0. 1  or  less. 

These  curves  serve  as  a  point  of  verification  that  the  numerical  approximations  of 
the  finite  element  technique  are  solving  the  governing  equations  correctly.  In  addition,  a 
point  of  validation,  indicating  that  the  correct  equations  are  being  solved,  is  desired.  To 
that  end,  the  model  was  exercised  over  several  specific  frequency-depth  combinations 
corresponding  to  recent  experimental  work  performed  at  the  Naval  Research 
Laboratory  [32],  Figure  27  reproduced  from  [32],  shows  a  schematic  of  the  experimental 
setup. 


48 


Arb,  Waveform 


Microphone  1 


air 

iS  cm 

Microphone  2 

water 

source 

;\ 

1  \  20  cm 

36  degi^ 

t  

3  mm 

\ 

hydrophone 


Figure  27.  Schematic  Diagram  of  the  Experimental  Setup  Used  in  the  NRL 

Experiments 

Source  [32]:  D.  C.  Calvo,  M.  Nicholas  and  G.  J.  Orris,  “Experimental  verification  of 
enhanced  sound  transmission  from  water  to  air  at  low  frequencies,”  Journal  of  the 
Acoustical  Society  of  America,  vol.  134,  no.  5,  pp.  3403-3408,  Nov.  2013. 


Two  airborne  microphones  were  placed  at  locations  on  and  off  axis  above  a  small 
spherical  transducer  located  a  small  depth  below  the  air-water  interface  of  an  acoustic  test 
tank.  A  hydrophone  was  also  located  at  a  fixed  depth  relative  to  the  source,  and  the 
relative  magnitudes  of  the  measured  air  and  water  borne  pressures  were  recorded  as 
functions  of  frequency.  The  ratios  of  these  experimentally  measured  pressures  and  their 
theoretical  predictions  are  compared  to  the  values  predicted  by  the  COMSOL  model  in 
Figures  28  and  29.  The  figures  suggest  excellent  agreement  with  the  theoretical 
predictions.  The  NRL  experimental  data  is  an  important  verification  of  the  theory. 
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Figure  28.  Comparison  of  the  Ratios  of  Acoustic  Pressure 

In  air  to  acoustic  pressure  in  water  plotted  against  non-dimensional  source  depth. 
COMSOL  (Left)  and  NRL  data  (right  dots)  along  with  theoretical  predictions  (right 
lines). 
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Figure  29.  Comparison  of  the  Ratios  of  Acoustic  Pressures 

Measured  off-axis  to  those  measured  on-axis  plotted  against  non-dimensional  source 
depth.  COMSOL  (Left)  and  NRL  experimental  data  (right  dots)  along  with  theoretical 
predictions  (right  lines). 
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III.  ANALYSIS 


Examining  the  plots  comparing  the  Finite  Element,  Normal  Modes,  and  Parabolic 
Equation  methods  reveals  several  interesting  differences.  First,  the  normal  modes  method 
consistently  predicts  higher  levels  in  all  fields.  Second,  COMSOL  consistently  predicts 
more  interference  patterns  and  features  in  the  acoustic  field  than  do  the  other  two 
methods.  Third,  the  PE  method  shows  better  agreement  with  the  other  two  modeling 
approaches  as  both  the  depth  of  the  water  column  and  frequencies  of  study  are  increased. 
Finally,  the  trade-off  between  feature  prediction  and  computational  problem  size  is 
discussed. 

In  all  cases,  regardless  of  the  treatment  of  the  bottom  interface,  the  normal  modes 
based  solution  predicts  acoustic  pressure  levels  higher  than  the  other  two  methods. 
Although  there  are  no  explicit  material  loss  mechanisms  present  in  any  of  the  models, 
there  are  several  sources  of  potential  attenuation  in  both  the  finite  element  and  parabolic 
equation  approach  that  could  account  for  the  discrepancy  in  predicted  pressure  levels.  In 
the  finite  element  approach,  it  is  well  known  that  the  density  of  the  mesh  will  affect  the 
solution.  While  the  results  presented  in  this  work  took  great  care  to  ensure  adequate  mesh 
resolution  to  predict  an  accurate  field,  prior  efforts  have  indicated  that  levels  may  slightly 
increase  in  magnitude  even  after  the  heuristically  determined  ten  linear  elements  per 
wavelength  is  reached.  The  meshes  used  in  this  study  were  refined  to  twenty  quadratic 
elements  per  wavelength,  and  in  the  case  of  the  rigid  bottom,  there  is  excellent  agreement 
in  levels  between  the  COMSOL  results  and  the  normal  mode  results.  The  very  small 
disparity  in  magnitudes  between  the  two  methods  may  be  brought  into  better  agreement 
with  further  mesh  refinement.  However,  it  was  noted  that  there  were  no  significant 
differences  in  the  features  predicted  by  COMSOL,  and  the  sub  deciBel  differences  were 
not  determined  to  contribute  to  the  accuracy  of  the  solution.  Therefore,  this  is  merely  a 
hypothesis,  but  is  consistent  with  prior  results. 

The  acoustic  fields  predicted  by  the  finite  element  method  unfailingly  predict 
more  structure  and  “fringing”  than  the  other  two  methods.  In  the  case  of  the  rigid  bottom, 

there  are  only  slight  discrepancies  between  the  COMSOL  and  normal  modes  output.  In 
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this  case,  the  agreement  is  quite  good.  It  is  worth  noting  that  the  normal  modes  method 
does  not  actually  propagate  a  solution  over  range,  but  rather  seek  solutions  to  an 
eigenvalue  problem  with  the  appropriate  boundary  conditions  already  enforced.  This  is 
closer  to  the  finite  element  based  matrix  methods  than  is  the  parabolic  equation  approach. 
The  results  of  the  parabolic  equation  solution  show  large  discrepancies,  failing  to  predict 
nearly  all  of  the  features  captured  by  the  other  two  methods  at  any  range.  However,  this  is 
a  well-understood  limitation  of  the  code,  which  is  not  designed  for  rigid  bottom 
environments. 

In  the  case  of  the  sediment  bottom,  there  is  better  agreement  between  the  three 
methods,  but  COMSOL  still  predicts  additional  interference  patterns.  The  discrepancy 
between  COMSOL  and  the  normal  modes  approach  is  potentially  due  to  the  presence  of 
evanescent  waves  and  so  called  “leaky”  modes.  The  initial  normal  modes  code  has  a  hard 
limit  to  examine  the  contribution  of  only  the  trapped  modes,  whereas  COMSOL  has  no 
such  restriction.  Upon  switching  to  the  coupled  modes  code  Couple07,  significantly 
better  agreement  between  all  three  methods  was  observed.  The  discrepancies  between 
COMSOL  and  MMPE  are  again  likely  due  to  the  treatment  of  the  bottom  interface.  There 
is  a  much  smaller  mismatch  in  the  impedances  of  the  two  domains,  however  the 
bathymetry  is  still  very  shallow  and  the  density  mixing  functions  employed  by  this 
version  of  MMPE  are  likely  contaminating  the  solution.  This  effect  will  be  present  at  all 
ranges.  At  short  ranges,  it  is  unlikely  that  the  two  methods  would  exhibit  good  agreement 
due  to  the  far-field  approximations  used  by  the  parabolic  equation  method  to  generate  the 
initial  field.  It  is  also  worth  noting  that  while  both  the  finite  element  method  and  normal 
modes  method  seek  solutions  to  Equation  1  (albeit  through  very  different  approaches), 
the  parabolic  method  seeks  solutions  to  a  modified  form  of  the  Helmholtz  equation 
defined  in  Equation  7. 

As  previously  stated,  the  discrepancies  between  the  methods  when  the  bottom 
sediment  is  included  are  most  likely  due  to  the  treatment  of  the  bottom  interface  density 
discontinuity  within  the  parabolic  equation  framework.  When  the  test  case  documented  in 
the  JCA  paper  was  examined,  there  was  significantly  better  agreement  between  MMPE 
and  COMSOL.  This  environment  is  twice  as  deep  (200  m)  and  analyzed  at  four  times  the 
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frequency  (200  Hz)  with  the  source  placed  at  half  the  depth  (vice  three  quarters  of  the 
depth)  of  the  water  column.  As  such,  the  ratio  of  excitation  wavelength  to  water  column 
depth  is  much  smaller  than  the  bulk  of  the  cases  presented  in  this  work. 

The  finite  element  approach  predicts  a  richer  feature  set  in  the  pressure  field  than 
do  the  normal  modes  and  parabolic  equation  methods.  However,  the  solution  time  and 
computational  problem  size  are  significantly  larger  when  using  COMSOL  as  compared  to 
the  other  codes.  Both  the  normal  modes  and  MMPE  executable  will  run  on  a  large 
desktop  PC  (2.7  GHz  8  core  Pentium  with  72  GB  of  RAM)  in  well  under  a  minute. 
Comparable  environments  require  several  minutes  to  solve  within  COMSOL.  It  is 
therefore  important  to  be  cognizant  of  the  accuracy  requirements  of  the  desired  numerical 
predictions.  At  higher  frequencies  and  deeper  bathymetries,  the  additional  features 
predicted  by  the  finite  element  method  may  be  unwarranted  due  to  the  computing  time 
necessary  to  generate  them.  Similarly,  if  the  overall  aim  of  study  is  to  predict  pressure 
based  features,  such  as  transmission  loss  alone,  the  accuracy  afforded  by  parabolic 
equations  or  normal  modes  is  likely  sufficient.  Here  however,  because  of  interest  in  the 
acoustic  intensity  field,  which  depends  on  both  the  pressure  field  and  its  gradient, 
additional  time  may  be  required  to  predict  the  necessary  features  accurately.  This  is  most 
notable  in  examination  of  Figure  18  where  the  general  shape  and  sign  of  the  intensity 
field  is  consistent  between  all  three  methods,  but  there  are  striations  predicted  only  by 
COMSOL. 

The  SLFT  processing  technique  allows  for  a  detailed  examination  of  the 
propagating  wave  types  for  arbitrary  bathymetries.  In  typical  experimental  data  collection 
events,  vertical  line  arrays  (VLA)  are  employed  in  order  to  gain  as  large  an  aperture  as 
possible.  VLAs  are  also  attractive  due  to  their  relative  ease  of  deployment,  naturally 
working  in  concert  with  Earth’s  gravitational  field.  Horizontal  arrays  resting  along  the 
ocean  seabed  are  less  common,  but  are  the  closest  physical  analog  to  the  digital  technique 
demonstrated  herein.  In  either  case,  a  long  array  is  desired  in  order  to  provide  adequate 
wavenumber  separation  during  data  processing.  For  relatively  shallow  water  littoral 
environs,  the  depth  of  the  water  column  may  preclude  sufficiently  long  aperture  sizes  to 
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examine  low  frequency,  large  wavenumber  propagation.  Employment  of  numerical 
models  and  synthetic  aperture  processing  therefore  offers  several  advantages. 

The  results  of  explicitly  modeling  the  air-water  interface  have  confirmed  recently 
published  theories  regarding  the  transparency  of  said  interface  [32],  and  conditions  under 
which  classical  theorems  hold.  In  most  practical  applications,  the  product  of  the  source 
depth  and  the  wavenumber  generated,  kz,  is  of  order  unity,  and  classical  theories 
describing  the  reflection  and  transmission  coefficients  are  adequate  to  describe  both 
propagating  and  evanescent  waves  [23].  However,  in  a  limited  number  of  cases,  such  as 
infrasonic,  or  very  shallow  source  depths,  the  product  of  kz  may  drop  to  a  value  of  order 
0.1  or  less.  Under  these  circumstances,  the  mathematical  treatment  of  Godin  [33]  is 
necessary  to  capture  all  the  governing  physics.  COMSOL  was  used  to  verify  this  theory 
and  validate  against  measure  data. 
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IV.  CONCLUSIONS 


Three  techniques  have  been  applied  to  the  problem  of  acoustic  propagation  in  a 
shallow  water  waveguide  at  low  frequency.  For  the  canonical  case  of  an  ideal  waveguide 
(pressure  release  and  rigid  boundaries),  there  is  excellent  agreement  between  the  theory 
of  normal  modes  and  the  finite  element  methods  predicted  using  the  commercial  software 
package  COMSOL.  There  are  significant  differences  in  the  predictions  of  the  parabolic 
equation  based  code  MMPE  that  may  be  attributed  to  the  implementation  of  the  blending 
functions  requisite  for  bottom  treatment.  With  the  inclusion  of  a  non-rigid,  sediment 
bottom,  idealized  as  a  second  acoustic  media  without  shear  or  attenuation,  there  is  much 
better  agreement  between  all  three  methods,  though  there  are  additional  features 
predicted  by  COMSOL.  This  is  attributed  to  the  method  by  which  COMSOL  treats  the 
scattering  interactions  from  boundaries  and  the  ability  of  the  finite  element  method  to 
capture  evanescent  waves.  The  prediction  of  these  features  comes  at  a  significant 
computational  cost  relative  to  other  methods,  but  is  necessary  to  predict  striations  in  the 
intensity  field  at  short  ranges.  Further,  the  intensity  features  present  in  the  second 
acoustic  media  are  explicitly  modeled.  In  cases  where  sensors  lie  at  or  below  the  seabed, 
the  accuracy  of  the  finite  element  method  warrants  the  additional  run-time.  In  many 
instances,  these  additional  near-field,  or  sub  seabed  features  in  the  intensity  fields 
predicted  by  the  finite  element  method,  may  not  be  of  great  importance.  Under  these 
circumstances,  the  PE  method  provides  adequate  resolution  and  allows  for  rapid 
parametric  model  evaluation.  In  cases  where  long  range  intensity  features  are  desired,  the 
well-known  phasing  issue  [15]  in  the  current  MMPE  code  may  present  sufficient  issues 
with  feature  prediction  that  a  normal  modes  code,  such  as  Couple07  is  recommended. 
However,  recent  developments  with  the  MMPE  model  have  begun  to  mitigate  this  issue 
to  improve  its  accuracy  [34].  Furthermore,  uncertainty  in  the  knowledge  of  the  ocean 
environments  at  long  range  may  overwhelm  such  errors  in  practice.  At  exceedingly  low 
frequencies,  a  full  treatment  of  the  air-water  interface  is  necessary  to  capture  the 
transmission  of  acoustic  energy  across  this  high-contrast  boundary. 
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